% This runs ResistancABMFinal.m multiple times, sweeping desired variables

%% Set Up
% Sweep Settings
OnlyGoodData = 0; % Set to 1 if you want cases filter to those deemed "Good"
PlotRaw = 0; % Set to 1 to plot raw data for each Good case
LogisticOrNot = 1; % Set to 1 for Logistic Regression, 0 for Rolling Average
LogOrNot = 1; % Set to 1 for logged peak participation, 0 for not
MatchUpTo = .5; % When calculating error, only calculate up to this value for X
SizeErrorGoal = .8; % Maximum error for peak participation when determining best case
NumberOfRuns = 500; % Number of runs for each parameter case
CountTimedOutAsLoss = 1; % Set to 1 to count timed out runs as losses

%Get NAVCO 1.2 Data
HistoricalData = GetNAVCOData(4,LogisticOrNot,LogOrNot,0,[],0,0); %Nonviolent Data, Logistic, Log, No plotting
PopulationBins = [-1:.5:1.5]; % Bins for Peak Participation Size (Only works with Log Option)
temp = GetRollingAverage(HistoricalData.XRaw,HistoricalData.YRaw,PopulationBins,1,.25,1);
HistoricalData.BinNumber = temp.nsize;clear temp; % Gets sum of peak participation cases for each bin

%% Declare Simulation Parameters
% Simulation Variables
sim.PauseTime = 0; % Seconds to pause each time step
sim.PlotOn = 0; % 1 to make plots, 0 for no plots
% Model Parameters
clear p
p.MaxSteps = 200; % Maximum number of time steps
p.LatticeX = 40; % Number of spaces in graph in X direction
p.LatticeY = 40; % Number of spaces in graph in Y direction
p.torus = 1; % Should lattice connect at sides? 1 Yes, 0 No
p.DelayStartMax = 5; % Maximum number of time steps an agent is delayed from acting intially
p.ReorderAgentsParam = 1; % Set to 1 if order of agents is randomly reordered
p.PercentFillCitizens = 70; % Percent Lattice is filled with agent type
p.PercentFillActivists = nan; % Percent Lattice is filled with agent type
p.PercentFillRev = 0; % Percent Lattice is filled with agent type
p.PercentFillPolice = 4; % Percent Lattice is filled with agent type
p.PercentFillPillars = .85; % Percent Lattice is filled with agent type
p.vision = 4; % Number of spaces an agent can see in each direction
p.MaxJailTerm = 10; % Max Jail Term in time steps
p.StartingGovernmentLegitimacy = .56; % Starting value (between 0 and 1) for government legitimacy
p.ChanceFindNVResistor = 40;
p.ChanceTargetNonviolent = 25; % Percent chance police will target a nonviolent agent
p.ChanceKillNonviolent = 10; % Percent chance police will attempt to kill nonviolent agent
p.BackfireCoefficient = .99; % Factor by which government legitimacy decreases
p.f = .0706; % Threshold used for determinging behavior
p.ProtestCycle = 7; % Frequency at which a nonviolent agent can protes
p.ProtestDuration = 1; %Number of time steps a nonviolent agent can protest continuously
p.nNV = 1; % Number of nonviolent people needed to persuade agent to protest
p.PeerPressureNumber = 3.3884; % # of Protestors for PeerPressure to equal .5 (and vision = 7)
p.PercentFickle = [100 75 50 25 0]; % Percent of Citizens who join and rejoin resistance based on seeing protest. The nonfickle stick remain regardless;
p.PercentImmediateProtest = 0; % %Percent of Citizens who will immediately protest after joining
p.DefectThreshold = .0547; % Mean defect threshold assigned to pillars and police
p.DefectThresholdStdDev = 0; % Defect Threshold does not vary pillar to pillar
p.NonviolentSuccessPercent = 1; % Percent of pillars required to defect as victory condition
p.PillarProxStrategy = 0; % Set to 0, 1, or 2 for different activist strategies
p.ActivistSearchVision = 10; % Vision an activist can look for pillars if Pillar Prox Strategy is 2
% Model Variation Parameters (i.e. parameters for variations run to run)
PopulationMean = .8; % Mean value for PercentFillActivists
PopulationStdDev = .3; % Standard Deviation value for PercentFillActivists
p.R2RDefectThresholdMin = .03; % Minimum value for Defect Threshold
p.R2RDefectThresholdSTD = .07; % Value multiplied by random number squared for Defect Threshold
p.R2RDefectThresholdMax = .3; % Maximum value for Defect Threshold
p.R2RNVSuccessPercentMin = 1; % Minimum value for NV Success Percent
p.R2RNVSuccessPercentSTD = 16; % Value multiplied by random number squared for NV Success Percent
p.R2RNVSuccessPercentMax = 80; % Maximum value for NV Success Percent
p.R2RPillarProxStrategyPercent = 0; % Value used to determine each run's PillarProxStrategy 

% Create Matrix of All Possible Variables
[m, ncases] = param2mat(p);

%% Runs Cases

%Sets Values to Be Varied each run of a given case
vary = randn(NumberOfRuns,1)*PopulationStdDev+PopulationMean; %Varies initial activist density
vary(vary<=.1) = .1;
vary2 = (randn(NumberOfRuns,1)).^2; % Needed For Chi Squared Distribution for NV Success Percent
vary3 = (randn(NumberOfRuns,1)).^2; % Needed For Chi Squared Distribution for Defect Threshold
vary4 = rand(NumberOfRuns,1); %Needed for Uniform Distribution for Pillar Proximity Strategy
    
% Initialization
tic % starts clock for recording how long the simulations take
clear out paramtemp LargestProtest wins ts LargestProtest2

% Run Simulations
for ii =  1:ncases
    paramtemp = m2instance(m,ii); % Selects Case from variable M

    % Creates Run to Run Variations
    vary2val = vary2*paramtemp.R2RNVSuccessPercentSTD + paramtemp.R2RNVSuccessPercentMin; % Chi-squared distribution
    vary2val(vary2val>p.R2RNVSuccessPercentMax)= p.R2RNVSuccessPercentMax; % Sets max value
    vary3val = vary3*paramtemp.R2RDefectThresholdSTD + paramtemp.R2RDefectThresholdMin; % Chi-squared distribution
    vary3val(vary3val>p.R2RDefectThresholdMax)= p.R2RDefectThresholdMax; % Sets max value
    vary4val = vary4 < paramtemp.R2RPillarProxStrategyPercent; % Set to 1 or 0
    
    % Runs Each Case Multiple Times
    %parfor jj = 1:length(vary) % Uses Parallel Computing Toolbox for speed
    for jj = 1:length(vary) % Does not require Parallel Computing Toolbox
    
        % Creates Final Parameter set
        temp = paramtemp;
        % **** The can be commented out to use non-varied parameter, i.e.
        % the same value for all runs ***
        temp.PercentFillActivists = vary(jj);
        temp.NonviolentSuccessPercent = vary2val(jj);
        temp.DefectThreshold = vary3val(jj);
        %temp.PillarProxStrategy = vary4val(jj);
        
        % Runs Model
        out = ResistanceABMFinal(sim,temp);
        
        % Gets Needed Results from Each Model
        LargestProtest(ii,jj)=out.MaxInProtestPercentage; % Saves Largest Peak Protest Size
        wins(ii,jj)=out.victory; % Saves win condition
        %ts(ii,jj) = out.nTimeSteps; % Saves number of steps when simulation finishes
    
    end
    casecounter = ii % Output to the screen when case complete
end
toc % ends clock
save TempSweep.mat % saves data

%% Post Process

% Determine Good Data
clear GoodData
if OnlyGoodData
    % Determines if there is more than one step from losses to wins
    for ii = 1:ncases
        [blah ind] = sort(LargestProtest(ii,:),2);
        tempwins = wins(ii,ind);
        GoodData.switches(ii) = sum(abs(diff(tempwins)))>1;
    end
    GoodData.switches = GoodData.switches';
    GoodData.total = GoodData.switches; % Other metrics could be added
else
    GoodData.total = true(1,ncases); % Sets all cases as Good Data
end
ind2plot = find(GoodData.total ==1); %Find Indicies of Good Data

% Make Plots of wins, losses, and timed out runs
if PlotRaw
    for ii = 1: length(nonzeros(GoodData.total))
        % Sorts data to make easier to read plot
        [blah ind] = sort(LargestProtest(ind2plot(ii),:),2);
        tempwins = wins(ind2plot(ii),:);
        % Plots
        figure
        plot(blah,tempwins(ind),'-x')
    end
end

% Gets Probability curves using Logistic Regression or Rolling Average
% Method
% *** Important Note ***, While the paper only comments on Logistic
% Regressions, another option is to create a plot based an average or
% rolling average to determine probability of success for a given protest
% size.
% Initialize variables
ProtestVector = [];
WinsVector = [];
CIlow = [];
CIhigh = [];
plusminus = .25; % Sets the x-value plus and minus around each x point to be averaged
nsize = [];
for ii = 1:ncases
    % Determine total number of runs
    CompletedRuns = wins(ii,:)>(0-CountTimedOutAsLoss); % Depending on CountTimedOutAsLoss, determines which runs to include in probability of success
    if sum(CompletedRuns)>0
        if LogisticOrNot % Use Logistic Regression Method
            xpoints = [-1:.1:1.5]; %Sets X-values to get datapoints (Must change for nonlogged case)
            temp = GetLogisticRegression(LargestProtest(ii,CompletedRuns)',[wins(ii,CompletedRuns)==2]',xpoints,LogOrNot,0);
        else % Use Rolling Average Method
            xpoints = [-1:.5:1.5]; %Sets X-values to get datapoints (Must change for nonlogged case) 
            temp = GetRollingAverage(LargestProtest(ii,CompletedRuns),wins(ii,CompletedRuns)==2,xpoints,LogOrNot,plusminus,0);
        end
    else
        % If no good runs, then it sets all values to "nan", i.e. not a
        % number
        xpoints = [-1:.1:1.5];
        temp.Y = xpoints.*nan;
        temp.CIlow = temp.Y;
        temp.CIhigh = temp.Y;
    end
    % Uses Average method to get the number of runs for specified peak
    % participation bins 
    temp2 = GetRollingAverage(LargestProtest(ii,CompletedRuns),wins(ii,CompletedRuns)==2,PopulationBins,LogOrNot,.25,LogOrNot);
    % Places temporary values into array with data from all runs
    WinsVector(ii,:) = temp.Y;
    CIlow(ii,:) = temp.CIlow;
    CIhigh(ii,:) = temp.CIhigh;
    ProtestVector(ii,:) = xpoints;
    nsize(ii,:) = temp2.nsize;
end

% Calculates Probabilty of Success and Outputs to Command Window
ProbSuccess = sum(wins==2,2)./sum(wins>(0-CountTimedOutAsLoss),2)*100

% ** Plotting All Data (Modify labels, title, legend as needed) **
% Make Population Graph (Unused in paper)
figure
% Normalizes number of runs for historical and model data
HistoricalData.PopScaled = HistoricalData.BinNumber./sum(HistoricalData.BinNumber);
PopScaled = nsize./repmat(sum(nsize,2),1,size(nsize,2));
% Plots
plot(PopulationBins,HistoricalData.PopScaled,'-ok',PopulationBins',PopScaled(GoodData.total,:)','-x')
ylim([0 .5])
xlim([-1 1.5])
legend('location','NorthEast')
xlabel('Peak Participation log10(%)') 
ylabel('Number of Cases')
% Make Success Graph 
figure
PlotWithCI = 1;
if PlotWithCI
    % Plot with Confidence Intervals
    handles = plot(HistoricalData.X,HistoricalData.Y,'-k',HistoricalData.X,HistoricalData.CIlow,'-.k',HistoricalData.X,HistoricalData.CIhigh,'-.k',ProtestVector(GoodData.total,:)',WinsVector(GoodData.total,:)','-',ProtestVector(GoodData.total,:)',CIlow(GoodData.total,:)','-.',ProtestVector(GoodData.total,:)',CIhigh(GoodData.total,:)','-.');
    % Modifies colors so all somewhat dark
    ColorList = {[0 0 1],[1 0 0],[0 .5 0],'m',[.6 .3 0],[.5 .5 .5],'y'};%Blue, Red, Dark Green, Magenta, Dark Orange, Gray, Yellow
    for ii = 1:ncases
        set(handles(ii+3),'color',ColorList{ii})
    end
    % Modifies Matching Colors for CI
    for ii = 1:length(WinsVector(:,1))
        SolidInd = 3+ii;
        %CIInd = 3 + length(WinsVector(:,1)) + 1 + (ii-1)*2;
        set(handles(SolidInd+length(WinsVector(:,1))),'Color',get(handles(SolidInd),'Color'))
        set(handles(SolidInd+length(WinsVector(:,1))*2),'Color',get(handles(SolidInd),'Color'))
    end
else
    % Plots without Confidence Intervals
    handles = plot(HistoricalData.X,HistoricalData.Y,'-k',ProtestVector(GoodData.total,:)',WinsVector(GoodData.total,:)','-');
    % Modifies colors so all somewhat dark
    ColorList = {[0 0 1],[1 0 0],[0 .5 0],'m',[.6 .3 0],[.5 .5 .5],'y'};%Blue, Red, Dark Green, Magenta, Dark Orange, Gray, Yellow
    for ii = 2:length(handles)
        set(handles(ii),'color',ColorList{ii-1})
    end
end
ylim([0 100])
xlim([-1 1.5]) % For Log Case
xlabel('Peak Participation log10(%)')
ylabel('Probability of Success')
legend('location','southeast')
% Select Correct Legend/Title As Needed
%legend off
%legend('Historical Data', 'Success Percent = 1%','Success Percent = 25%','Success Percent = 50%','Success Percent = 75%','Success Percent = 99%','location','Northwest')
%legend('Historical Data', 'Defect Threshold = .025', 'Defect Threshold = .05', 'Defect Threshold = .1', 'Defect Threshold = .2', 'Defect Threshold = .3', 'Defect Threshold = .4', 'Defect Threshold = .5')
%legend([handles(1) handles(2) handles(4) handles(7) handles(5) handles(8) handles(6) handles(9)],'Historical Data', 'Historical CI','Model Run 1','Model Run 1 CI','Model Run 2','Model Run 2 CI','Model Run 3','Model Run 3 CI')
%legend('Historical Data','Max Steps = 150','Max Steps = 200','Max Steps = 250','Max Steps = 300')
%title('Ignoring Timed Out Cases')
%legend('Historical Data','Pillar Prox Strategy = 0', 'Pillar Prox Strategy = 1')
%legend('Historical Data','Pillar Prox Strategy = 0', 'Pillar Prox Strategy = 1','Pillar Prox Strategy = 2')
%legend('Historical Data','Pillar Prox Strategy = 0', 'Pillar Prox Strategy = 1','Pillar Prox Strategy = 2')
%legend('Historical Data','Historical Data 95% Conf. Interval','','Pillar Prox Strategy = 0','Pillar Prox Strategy = 1','Pillar Prox Strategy = 2','Pillar Prox Strategy = 0 95% Conf. Interval','Pillar Prox Strategy = 1 95% Conf. Interval','Pillar Prox Strategy = 2 95% Conf. Interval')
%h = get(gca,'Children');
%set(gca,'Children',[h(1) h(2) h(3) h(4) h(7) h(5) h(8) h(6) h(9) h(10) h(11) h(12)])
%legend('Historical Data','% Immediate Protest = 0%','% Immediate Protest = 25%','% Immediate Protest = 50%','% Immediate Protest = 75%','% Immediate Protest = 100%')
%legend('Historical Data','% Fickle = 100%','% Fickle = 75%','% Fickle = 50%','% Fickle = 25%','% Fickle = 0%')
%legend('Historical Data','% Committed = 0%','% Committed = 25%','% Committed = 50%','% Committed = 75%','% Committed = 100%')

% Find Best Case for Matching Historical Data
MatchUpToVector= xpoints<=MatchUpTo; %As diffferent cases have different peak sizes, only match date up to a certain portion of the fit
InterpHistoricalWins = interp1(HistoricalData.X,HistoricalData.Y,xpoints,'pchip','extrap'); % Interpolates Historical data to have Y values for the same X values as model data (so they can be compared)
% Calculates Normalized Size Error (and outputs to Command Window
SizeErrorRaw = sum((PopScaled-repmat(HistoricalData.PopScaled,ncases,1)).^2,2);
SStotPop = sum((HistoricalData.PopScaled-mean(HistoricalData.PopScaled)).^2);
SizeError = SizeErrorRaw/SStotPop
% Calculates Normalized Success Error (and outputs to Command Window
SuccessErrorRaw = sum((WinsVector(:,MatchUpToVector)-repmat(InterpHistoricalWins(MatchUpToVector),ncases,1)).^2,2);
SStot = sum((InterpHistoricalWins(MatchUpToVector)-mean(InterpHistoricalWins(MatchUpToVector))).^2);
SuccessError = SuccessErrorRaw/SStot

% Determine Array Index for Best Case
if sum(SizeError<=SizeErrorGoal)>0
    % If any cases have size error less than goal, then pick the case with
    % lowest success error of these cases
    ErrorTemp = SuccessError;
    ErrorTemp((GoodData.total==0)' | [(SizeError>SizeErrorGoal)])=nan;
    [MinError MinIndex] = min(ErrorTemp);
else
    % If no cases have size error less than goal, pick the case that has
    % lowest size error
    [MinError MinIndexISH] = min(SizeError(GoodData.total,:));
    MinIndex = ind2plot(MinIndexISH);
end

% Outputs Information about Best Matching Case to MATLAB Command Window
BestParam = m2instance(m,MinIndex)
BestSuccessError = SuccessError(MinIndex)
BestSizeError = SizeError(MinIndex)
BestMatchProbSuccess = ProbSuccess(MinIndex)

% Plots Best Case: Participation Size Histogram
figure
hh = bar(PopulationBins,[HistoricalData.PopScaled', PopScaled(MinIndex,:)'],'k');
set(hh(2),'FaceColor',[0 0 1],'EdgeColor',[0 0 1]); % Adjusts color of model data bar
ylim([0 .4])
xlim([-1.25 1.75]) % For Logged X-axis
%xlim([0 3]) % For Linear X-axis
legend('Historical Data','Model Data','location','NorthEast')
xlabel('Peak Participation log10(%)') % Need to change if not log plot
ylabel('Normalized Number of Cases')

% Plots Best Case: Probability of Success for Given Participation Size
figure
hlr = plot(HistoricalData.X,HistoricalData.Y,'-k',HistoricalData.X,HistoricalData.CIlow,'-.k',HistoricalData.X,HistoricalData.CIhigh,'-.k',ProtestVector(MinIndex,:)',WinsVector(MinIndex,:)','-b',ProtestVector(MinIndex,:)',CIlow(MinIndex,:)','-.b',ProtestVector(MinIndex,:)',CIhigh(MinIndex,:)','-.b');
ylim([0 100])
xlim([-1 1.5]) % For Logged X-axis
%xlim([0 3]) % For Linear X-axis
legend([hlr(1) hlr(2) hlr(4) hlr(5)],'Historical Data','Historical CI','Model Data','Model CI')
legend('location','southeast')
xlabel('Peak Participation log10(%)') % Need to change if not log plot
ylabel('Probability of Success')


